ILF <- function(X,W = rep(1/NROW(X),NROW(X)),ot,n,R = NULL){
  # This function calculates the Inverse Lorenz Function and related objects.
  # Output of this function is intended to be used to sketch contours of the ILF,
  # which are the alpha-Lorenz curves.
  #
  # Output: lmap       = Lorenz map estimates at each vector rank sampled from U([0,1]^2)
  #         estimate   = The ILF on a grid, i.e., the empirical CDF of lmap
  #         Gini       = bivariate Gini calculated as plug-in using lmap
  #         shares     = result of calculating lmap for specified ranks R
  #         quantile   = power diagram representing the vector quantile of data
  #         constraint = max difference between areas of the cells and constraint
  #                      to roughly test if constraints are satisfied
  #         flag       = counts how many empty cells resulted
  #         time       = time elapsed
  #
  # Inputs: X  = bivariate allocation (data)
  #         W  = sampling weights of X, default is equal weights (1/NROW(X))
  #         ot = pwd output of "vquantile.R", includes power diagram weights
  #         n  = size of sample to take from the uniform distribution on [0,1]x[0,1]
  #         R  = (optional) k-by-2 matrix of vector ranks to produce desired vector shares
  #               This is necessary to generate figure 6 in the paper            
  #
  # Summary: Step 1: Calculate L(R_i) for each i = 1,...,n and R_i from U[0,1]^2
  #                  (Done in parallel)
  #          Step 2: Calculate the ecdf of L(R_i). 
  #
  
  # Step 0: Preliminaries
  start <- Sys.time()
  df <- data.frame(X,W)
  colnames(df) <- c("X1","X2","weight") 
  df <- aggregate(weight ~ X1 + X2, df, sum) # aggregate based on duplicates
  W <- df$weight/sum(df$weight)
  mu <- cbind(df[,1]%*%W, df[,2]%*%W)
  X <- cbind(df[,1]/mu[1],df[,2]/mu[2])
  pwd <- ot$pwd
  
  # Step 1
  N <- length(W)
  cl <- makeCluster(detectCores()-1, type = "SOCK")
  clusterExport(cl=cl,varlist = c('Lmap','X','W','pwd','st_area','st_intersection','st_sfc','st_polygon','rowProds','N'),envir=environment())
  # Refer to "Lmap.R" for function details
  LR <- do.call(rbind,parLapply(cl,1:n, function(i){Lmap(cbind(runif(1),runif(1)),X,W,pwd$cells)}))
  stopCluster(cl)
  
  LR <- LR[complete.cases(LR),]
  
  # Step 2: Calculate the ecdf of L(R_i).
  # ecdf of LR is calculated on a fine grid.
  grid <- as.matrix(CJ(x = seq(0,1,0.005), y = seq(0,1,0.005)))
  Z <- apply(grid,1,function(r){eCDF(LR,r[1],r[2])})
  Z <- matrix(Z, nrow = length(seq(0,1,0.005)))
  
  # Calculate bivariate Gini coefficient by plug-in using L(R_i) sample
  G <- 1 - 2*(sum(colMeans(LR)))
  
  # Calculate vector cumulative shares for the bottom vector ranks specified in R
  if (is.null(R)){shares <- NULL}else{
  shares <- sapply(1:NROW(R), function(i){Lmap(c(R[i,]),X,W,pwd$cells)})
  }
  
  end <- Sys.time()
  list_return <- list(LR,Z,G,shares,pwd,ot$constraint,ot$empty,end - start)
  names(list_return) <- c('lmap','estimate','Gini','shares','quantile','constraint','empty','time')
  return(list_return)
}